function hfile = vmr_InhomogeneityCorrect(hfile, opts)
% VMR::InhomogeneityCorrect  - attempt automatic inhomogeneity correction
%
% FORMAT:       [vmr = ] vmr.InhomogeneityCorrect([opts])
%
% Input fields:
%
%       opts        optional struct with settings
%        .mask      either 3D uint8/logical data or VMR object with preseg
%                   if omitted, try automatic mask detection
%        .model     either of 'log', {'mult'}
%        .order     polynomial order, default 3 (valid: 2 through 7)
%        .trf       if ACPC and TAL objects given, use apriori selection
%
% Note: if no mask is given and the trf field is set, this function uses
%       a prior file compiled out 20 WM segmentations (0.5mm coordinates)

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if nargin < 2 || ...
    numel(opts) ~= 1 || ...
   ~isstruct(opts)
    opts = struct;
end
if ~isfield(opts, 'mask') || ...
    isempty(opts.mask)
    opts.mask = [];
end
if numel(opts.mask) == 1 && ...
    isBVQXfile(opts.mask, 'vmr')
    mbc = bvqxfile_getcont(opts.mask.L);
    opts.mask = mbc.VMRData;
    if ~isequal(size(opts.mask), size(bc.VMRData))
        opts.mask = [];
    end
    if ~isa(opts.mask, 'uint8')
        opts.mask = uint8([]);
    end
    opts.mask(opts.mask < 226) = 0;
    opts.mask(opts.mask > 225) = opts.mask(opts.mask > 225) - 225;
    um = unique(opts.mask(:) + 1);
    ur = uint8(1:max(um));
    ur(uo) = 1:numel(um);
    opts.mask = ur(opts.mask);
end
if ~isfield(opts, 'model') || ...
   ~ischar(opts.model) || ...
   ~any(strcmpi(opts.model(:)', {'l', 'log', 'm', 'mult'}))
    opts.model = 'mult';
else
    opts.model = lower(opts.model(:)');
end
if ~isfield(opts, 'order') || ...
    numel(opts.order) ~= 1 || ...
   ~isa(opts.order, 'double') || ...
    isnan(opts.order) || ...
   ~any((2:7) == opts.order)
    opts.order = 2;
end
if ~isfield(opts, 'trf') || ...
   ~iscell(opts.trf) || ...
    numel(opts.trf) ~= 2 || ...
    numel(opts.trf{1}) ~= 1 || ...
   ~isBVQXfile(opts.trf{1}, 'trf') || ...
    numel(opts.trf{2}) ~= 1 || ...
   ~isBVQXfile(opts.trf{2}, 'tal')
    opts.trf = {};
else
    opts.trf = opts.trf(:)';
end

% and data size
vs = size(bc.VMRData);
vtx = vs(2) * vs(3);
vty = vs(1) * vs(3);
vtz = vs(1) * vs(2);

% voxel selection thresholding
xi = find(ihd.numx >= (opts.voxselt * vtx));
yi = find(ihd.numy >= (opts.voxselt * vty));
zi = find(ihd.numz >= (opts.voxselt * vtz));
if any([numel(xi), numel(yi), numel(zi)] < (opts.order + 2))
    error( ...
        'BVQXfile:InternalError', ...
        'Voxel selection thresholding failed.' ...
    );
end

% find centers
ofx = -find(cumsum(ihd.numx(xi(1):end)) >= sum((ihd.numx)) / 2);
ofy = -find(cumsum(ihd.numx(xi(1):end)) >= sum((ihd.numx)) / 2);
ofz = -find(cumsum(ihd.numx(xi(1):end)) >= sum((ihd.numx)) / 2);
ofx = ofx(1);
ofy = ofy(1);
ofz = ofz(1);

% make fit
mnx = mean(ihd.meanx(xi));
mny = mean(ihd.meany(yi));
mnz = mean(ihd.meanz(zi));
px = polyfit(xi + ofx, ihd.meanx(xi) - mnx, opts.order);
py = polyfit(yi + ofy, ihd.meanx(yi) - mny, opts.order);
pz = polyfit(zi + ofz, ihd.meanx(zi) - mnz, opts.order);
sx = (1:vs(1)) + ofx;
sy = (1:vs(2)) + ofy;
sz = (1:vs(3)) + ofz;

% what data
if isequal(size(bc.VMRData), size(bc.VMRData16))
    vd = single(bc.VMRData16(:, :, :));
    v16 = true;
else
    vd = single(bc.VMRData(:, :, :));
    v16 = false;
end

% create correction matrix
if opts.model(1) == 'a'
    cvd = ...
        repmat(reshape(-polyval(px, sx), [numel(sx), 1]), [1, vs(2), vs(3)]);
    cvd = cvd + ...
        repmat(reshape(-polyval(py, sy), [1, numel(sy)]), [vs(1), 1, vs(3)]);
    cvd = cvd + ...
        repmat(reshape(-polyval(pz, sz), [1, 1, numel(sz)]), [vs(1), vs(2)]);
    vd = vd + cvd;
    vd(vd < 0) = 0;
else
    pv = mnx ./ (mnx + polyval(px, sx));
    pv(pv < opts.limit(1)) = opts.limit(1);
    pv(pv > opts.limit(2)) = opts.limit(2);
    cvd = osz(vd) .* repmat(reshape(pv, [numel(sx), 1]), [1, vs(2), vs(3)]);
    pv = mny ./ (mny + polyval(py, sy));
    pv(pv < opts.limit(1)) = opts.limit(1);
    pv(pv > opts.limit(2)) = opts.limit(2);
    cvd = cvd .* repmat(reshape(pv, [1, numel(sy)]), [vs(1), 1, vs(3)]);
    pv = mnz ./ (mnz + polyval(pz, sz));
    pv(pv < opts.limit(1)) = opts.limit(1);
    pv(pv > opts.limit(2)) = opts.limit(2);
    cvd = cvd .* repmat(reshape(pv, [1, 1, numel(sz)]), [vs(1), vs(2)]);
    vd = vd .* cvd;
end

% what result
if v16
    vd(vd > 32767) = 32767;
    vd = uint16(vd);
    bc.VMRData16 = vd;
else
    vd(vd > 225) = 225;
    vd = uint8(vd);
    bc.VMRData = vd;
end

% set in output
bvqxfile_setcont(hfile.L, bc);
